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Abstract 

5»J , A method for generating random (7(1) variables with Boltzmann dis- 

tribution is presented. It is based on the rejection method with transfor- 
mation of variables. High efficiency is achieved for all range of tempara- 
tures or coupling parameters, which makes the present method especially 
suitable for parallel and pipeline vector processing machines. Results of 
pH [ computer runs are presented to illustrate the efficiency. An idea to find 

such algorithms is also presented, which may be applicable to other dis- 
tributions of interest in Monte Carlo simulations. 

• f-n - 

'^ ■ Subject classification: 65C10, 81E25, 82A68. 
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1 Introduction. 

Monte Carlo numerical integration method, or Monte Carlo simulation, has 
been widely used in the numerical study of quantum field thcroies with lat- 
tice formalism and statistical mechanics of spin systems. In performing Monte 
Carlo simulations, one must generate sequences of random numbers with given 
probability distributions. Each random number is used to 'update' a spin or 
a gauge variable. Probability distributions which appear in such calculations 



have parameter dependences. These parameters carry the information of tem- 
perature and other thermodynamic quantities, and also that of the neighboring 
spin states. 

For a fixed probabihty distribution there are algorithms (see, for example, 
[pi and references therein) which may be very efficient in generating random 
numbers. In Monte Carlo simulations, however, we have fluctuations in the 
neighboring spins, which results in changes of the parameters within a single 
program. One faces the problem of finding an algorithm suitable to a class of 
distributions parametrized by the parameters. These parameters are changed 
over a wide range, so that we must find a method which maintain high efhciency 
for all range of parameters, especially in the limit where a parameter tends to 
oo and the distribution becomes singular. 

Another aspect which we consider, is the efficiency in parallel or pipeline 
vector processing. As far as we know, much of the currently available vector 
processors work efficiently when there are no 'if-branches' in the program. 

We study a method based on a rejection method combined with a change 
of variables [Q , which is an approach that is widely used. In the case of rejec- 
tion method, to avoid 'if-branches' we have to fix the number of iterations of 
the rejection trials. It is particularly important to have high acceptance rate 
uniformly in the parameters. 

The aim of the present study is to find a suitable method for generating a 
sequence of random t/(l) numbers which we use to update site or link variables 
of a canonical ensemble for U{1) lattice gauge theories or f/(l) spin sytems, from 
the point of view mentioned above. An example of a program for our proposal 
is given in appendix O. 

In section we set up the problem. In section H we give our strategy for the 
solution, and in section we give an answer to the problem. We present the 
results of our efficiency tests of the algorithm in section ||. 

2 Random U{1) variable and rejection method. 

In the following we call a sequence of random U{1) numbers, a random L'^(l) 
variable: A random U{1) variable is a sequence of numbers (the angle variables) 

^1 , ^2 , 6*3 , • • ■ , (1) 

whose distribution 

Pi[0,e + d0])=fai0)d0 

is given by the density function 

fa{0) = NaeMacosie-0o)), 

where Na is a normalization constant. In the practical applications, the param- 
eter a is proportional to the inverse temperature 1/T or the inverse coupling 



l/g^, and both a and the constant 0o contain the effect of interactions with 
other sites or hnk variables. By the shift of variable 9' = 6 — do if a > , and 
9' = 9 — 9o — IT ii a < , we may assume without loss of generality that 9q = 
and a > . Hence 

fai9) = Naexp{acos9) , -n < 9 < tt , a > , 
1 



cx.p{a cos 9) d9. (2) 

The right hand side of eq. (||) is (modulo constant) an integral representation 
of modified Bessel function /q (a) . 

On computers we start with uniform random variables 

UJi , UJ2 , UJ3 , ■ ■ ■ , (3) 

with the probability distribution 

P{[uj, uj + duj]) = duj , < t^ < 1 . 

If we know an expression for a function X{u!) expressible as a computer program 
of small time consumption such that the sequence 

9,^X{lu,), i = l,2,3,---, 

is the random f/(l) variable (m), then there is in principle no problem. Such a 
function X{x) is formally given by 

X{x) = F-\x); 

F{t) = P{[-7r < 9 < t]) = f fa{9)d9. 

(Here and in the following, F^^ is the inverse function of F ; F~^{F{x)) — x .) 
Unfortunately, we do not know a suitable expression of X{x) for random U{1) 
variables. (Note that wc have the parameter a dependence.) 

Usually the rejection method is adopted to solve the problem. The rejection 
method, combined with transformation of variables, is defined as follows. Let 
f{9) be some approximate density function to the density function fa{9) . We 
assume that the density functions are continuous. Note that they are normalized 

r 

to satisfy / f{9)d9 — 1 . Suppose that there is a monotonic function h which 

satisfies 

/i(0) = -7r, /i(l) = 7r, 

and 

f{h{x))h'{x) = 1, 0<a;< 1, (4) 



where h' is the derivative of /i , h' — —— . (For the moment, we suppress possible 

ax 

parameter dependences of / and h .) Define a function g by 

Let cijj and a;' witli j = 1, 2, 3, • • • , be two sequences of independent uniform 
random variables as in (0). Define a subsequence 

uji =ujji, i = 1,2,3, ••• , 

of the sequence {ojj} by selecting the numbers j — ji that satisfy 

^j < gi^j) ■ 

Then the sequence 

h{oJi) , h{uj2) , /ilwa) , • • • , (7) 

is the random U{1) variable we are looking for. We call the rate of picking up 
LOi out of Uj the acceptance rate. To achieve high efficiency, the acceptance 
rate should be high. It can be shown that the acceptance rate is equal to R{a) 
in eq. (ph. The proof that the distribution of (M) is the one with the density 
function fa , and the proof that the acceptance rate is equal to R{a), are given 
in Appendix H. 

Note that the density functions are non-negative functions and that the 
integrated values of/ and fa are normalized to unity, hence < R{a) < 1 , which 
should hold for an acceptance rate. Note also that if / is a good approximation 
to fa , then the function g is almost flat. 

Our problem is to find a good / . 

3 Approximate distributions and the optimiza- 
tion of the acceptance rate. 

In order to keep high efficiency, we must choose / with high acceptance rate 
R{a). To illustrate the implications of this statement, let us first consider the 
simplest choice of / , the fiat distribution. Hereafter, we refer to this choice as 
the 'direct' method. 

The 'direct' method is defined by choosing / in cq. (0) to be a constant 
function; 

m = i- (8) 

ZTT 

h{x) = (2a;-l)7r. (9) 



The acceptance rate is 



R(a) = — — — . (10) 



For a near zero ('high temperature'), the acceptance rate is high; R{a) « 1 — a + 

3 

—a , a ^ 1 , while for large a ('low temperature'), the acceptance rate becomes 

very small; 

R{a) « , a > 1 . 

v27ra 

In the limit a ^ oo , the acceptance rate approaches zero. We have very slow 
effective generation of random U{1) variables for large values of a with the 
'direct' method. In other words, the improvement of the efficiency of generating 
random U{1) variables, which we measure by time lapse per a random variable 
generation, depends basically on the improvement of the acceptance rate. The 
reason that the acceptance rate is small for large a is that for large a the original 
distribution fa has a large peak at 6* = . The distribution becomes highly non- 
uniform and the flat (uniform) distribution / is not a good approximation. 

As noted in the section |l| (see also section ^, it is even more important 
to keep high acceptance rate uniformly in a when we use parallel processing 
machines. 

Our aim is to find a / which is a good approximation to the original distri- 
bution fa for all values of a , in particular, the a —* oo limit. In a sense, this is 
to find a family of distributions which interpolates the flat distribution (a = 0) 
and the delta function distribution [a — oo) expressible as a simple computer 
program. 

For large a, the density function fa has a sharp peak at 6* = . Therefore 
/ should be a very good approximation at — Q , while being not too bad an 
approximation a,t 6 — tt . Therefore it is desirable to have two free parameters 
for / . Also since fa is an even function, / is desired to be an even function. 
We also impose the condition that the corresponding function h in eq. (Q) has 
an analytic expression. The simplest choice satisfying these conditions is the 
following fa, 13 ; 

^"-^(^^- 2coshM)+2r ">"'^>-^' (") 

where Na^p is a normalization constant; 



Na,0 = < 



2ta.nh-\AB) 



"I' 



/?=!, (12) 



2tan-i(AB) 



where we put A — tanh — , and B = 

2 ' 



/3 + 1 



. The corresponding function h 



m eq. 



is; 



ha,p{x) = < 



- tanh" VS"4anh((2x - 1) tanh" VAB))) , (3 > 1 , 
a 

-tanh"^((2x- I) A) , (3=1, 

a 

2 

-tanh"^ (S"4an((2x- 1) tan"i(AB))) , 



-1</3<1, 

(13) 
where A and B are as above. 

The next step is to choose a = a{a) and /3 = (3{a) as functions of a . In 
principle, they should be chosen so as to optimize the acceptance rate R = R{a) . 
Here, we search for a solution that satisfies a condition that the minimum in the 
definition of R{a) (i.e. in the right hand side of eq. (^)) is achieved at 6 = . We 
impose this condition to avoid 'if-branches' in the resulting computer program. 
If the minimum is attained at different values oi 9 , we will need if-branches 
according to the different values of 9 . 

We have an argument that the optimal solution under this condition, which 
we shall refer to as the 'optimized cosh' method, is given by choosing a = a{a) 
and /3 = /3(a) in eq. dl3) to satisfy the following: 



1 



a(a) = V3a-1 , /3(a) = 2 - - , if a > a° . 

a 



cosh(7ra(a)) — 1 exp(2a) — 1 ^, , a(a)^ 
r-T2 = : P(a) = 



1 , if a° > a > a* 



(14) 



(15) 



Here a° and a* are positive constants satisfying a° > a*, uniquely determined 

by, 



exp(2a°) - 1 

a° 
exp(2a*) - 1 



cosh(7rV3a° - 1) - 1 



3a°-l 



Their numerical values are a 
eq. (H) is , for a > a* , 



Y' 

0.799 and 



(16) 
(17) 



i 5.04. The function g = ga in 

a{x) = exp{-aGa{ha{a),l3{a)ix))) (18) 



with 



Ga{e) = 1 - cos - i log ( 1 + — i— 

a \ 1 + /3(a) 



(cosh(a(a) 9) - 1) 



(19) 



For the parameter range of < a < a* , we have to take a limit a i with 
— -— fixed to 2 7r~^(exp(2a) — 1) . We have, in place of eq. (p^), 

h^{x) = -tan((2x- 1) tan"i(7r7)) , (20) 

7 

where 7 = 7(a) is 



7(a) = vr- Vexp(2a) - 1 , if < a < a* . (21) 

The function g = ga in eq. (||) is 

ga{x) ^exp{-aGa{h^(a){x))) (22) 

with 

G'„(6») = 1 -cos 6*- ilog(l + 7(a)26'2) . (23) 

The distribution is reduced to the Cauchy distribution: 

m = TT^^ (24) 

where 

^-y = ^TT 7- 

^ 2tan-i(7r7) 

See appendix H for the proof that these formulae correctly generate a random 
C/(l) variable, and arguments for our choice of the parameters. 

The acceptance rate R = R{a) for the 'optimized cosh' method is given by 

Ria) = ^-(-lf(-) , a>a\ 

^ ' 2Naexp{a){l + P{a)) ' " ' 

^^") a*>a>0, 



2 A^aexp(a) tan (7r7(a)) 



where Na and Na^p are defined in eq. (g) and eq. (Jl^), respectively. Note the 
high acceptance rate for both small a and large a: 

1 71 o 

R{a) « l--a+—a\ a<l. 



R(a) — > r=r- w 0.95, a — > 00 . 

21og(2 + V3) 

See Fig. 1 for the acceptance rate for full range of a. The acceptance rate for 
the 'optimized cosh' method keeps more than 90% for all values of a, including 
the 'zero temperature' limit a — > 00 . 



4 Proposed algorithm. 

The acceptance rate for the 'optimized cosh' method is high, but to obtain the 
parameter a{a) for a* < a < a" , one has to solve a transcendent equation 
eq. (^5|). Also, one has to use different formulae for < a < a* , a* < a < a° , 
and a" < a , which will cause 'if-branches' that will considerably lower the 
efficiency when using with vectorized processors. 

One may, for example, use Newton method to solve equations numerically, 
but here instead we approximate the function a{a) in eq. ( |l5|) directly by a 
function which is explicitly expressible on computer programs without using if- 
branches for all a, and designed to keep high acceptance rate. The if-branches 
are avoided in such a way that we have —1 < (3 < 1 for all a £ (0, oo). 

We shall give an example of realistic choice. 



1. Define a{a) and /3(a) by, 



a{a) = min{-\/a (2 — e) , maxjy^, 5{a)}} , 

a{af cosh(7ra(a)) - 1 

/3(a) = max{ , -— ^} - 1 , 

a exp(2a) — 1 



where 



(5(a) = 0.35 max{0 , a — a*} + 1.03 v^maxjO , a — a*} , 
and e = 0.001. a* = 0.798953686083986 is the constant defined by eq. (|l 
2. Define functions /ia,/3 and ga by 



ha,p{x) ^ - ta,nh M W -— — tan((2a:: - 1) tan ^(W (tanh— ))) 

and 

ga{x) = exp{-a Ga{ha(a),l3(a){x))) 

with 

Gje) = 1-C0S6I- -logf 1 + — -(cosh(a(a)6i) -1) 

a \ 1 + fj(a) 

3. Let ujj and lu', with j = 1, 2, 3, • • • , be two sequences of independent ran- 
dom variables uniformly distributing in [0, 1). Define a subsequence 

uji ^ ojji, i = 1,2,3, ••• , 

of the sequence {ujj} by selecting the numbers j = ji that satisfy 

^'j < 9a{uJj). 



The sequence 

is the random C/(l) variable, a sequence whose distribution is 

P{[e,9 + d9])^ Naexp{a cos e)d9, -tt < 9 < tt , a>0, 
1 



exp(acos6') d0 . 

We wiU refer to this choice as the 'proposed cosh' method. One can exphcitly 
check that the choice of parameters satisfies — 1 < P{a) < 1 for aU a g (0, oo). 
The functions /iq,,/3 and ga therefore are obtained from section |3|. One can 



see that this choice of parameters satisfies the conditions (32), (p3D, and (|34D, 
which impUes that the 'proposed cosh' method correctly generates C/(l) random 
variables. The acceptance rate R{a) for 'proposed cosh' method is given in 
Figure 1. Note that the acceptance rate is high for all values of a. The minimum 
of R{a) for 'proposed cosh' method is R{a = oo) « 0.88. 

5 Efficiency test. 

We will show our results of efficiency test for 'propsed cosh' method. 

In Monte Carlo simulations, we try to generate a number in the random 
[/(I) variable, for each trial update of a site or link variable. As explained in 
section 0, this is to generate two independent uniform random number ujj and 
Lu'j and decide whether to accept or reject h{ujj). We may alternatively decide to 
try at most n-times, namely to try uij, Wj+i, • • • , ujj+n-i, before deciding that 
this site variable was not updated. (We shall call this number n, the iteration 
number, and call this one set of trial, an update.) This will effectively improve 
the acceptance rate to 

i?„ = 1 - (1 - i?)" , (25) 

where R is the original acceptance rate. On the other hand, this will increase 
the time consumption per update by an average rate of roughly R~^ when 
n(l-i?)"< 1. 

In other words, one must consider the possibility of iterating an algorithm 
with low acceptance rate but high speed, such as iterating the 'direct' method. 
Improvement of the efficiency resides in a balance of high acceptance rate and 
simple (fast) program. 

Efficiency is a measure of average speed of producing random U{1) variable. 
From the above consideration on Monte Carlo updates, we may define the effi- 
ciency by the average time consumption for an update with iteration number n 
chosen so that the effective acceptance rate i?„ exceeds some fixed rate, say, 

Rn > 0.9 . (26) 



This definition is useful wfien considering Monte Carlo calculations. 

We note that this is not the only possible definition of efficiency. Alterna- 
tively, one may use the criteria to iterate until a variable is updated. In ap- 
pendix A of Q this criteria is adopted to measure the speed of their methods. 
However, except for the discreteness of iteration number (and modulo multipli- 
cation of a constant), the two criteria essentially measure the same quantity, so 
the conclusion will not change. 

The best possible choice depends on the machine to be used. We checked 
our choice in section |j by HITACH S820 in Computer Centre of University of 
Tokyo, which uses a pipelined vector processor. We measured the efficiency 
of generating random C/(l) variable defined as above, with average taken over 
4x 10^ updates. The 'proposed cosh' method (without iteration) has acceptance 
rate of more than 90% for a < 8 . We compared the speed of 'proposed cosh' 
method with the iterated 'direct' method for a < 8 , with the condition that the 
number of iteration is chosen to keep the acceptance rate to be more than 90%. 

The measured acceptance rates were in good agreement with the theoretical 
predictions in section || and in Figure 1 (within 0.1% acuracy). 

The results of the efficiency test is given in Figure 2, which shows that the 
'proposed cosh' method is efficient for all values of a, especially for large a. Note 
that for the 'direct' method, the CPU time increases as a increases, which is 
due to the decrease in the acceptance rate R{a) . Since R{a) — > as a ^ oo , the 
CPU time for the 'direct' method will tend to oo as a ^ cx3 . The acceptance 
rate for the 'proposed cosh' method is greater than 0.88647 for all values of a, 
consequently the CPU time is bounded for all a, and furthermore, the method 
is suitable for vector processors. 



6 Concluding remarks. 

We have shown a simple method for generating random C/(l) variables. Our 
method, based on rejection algorithm, achieves high acceptance rate and low 
time consumption for all values of the parameter a, including the 'low tempara- 
ture limit', a —> oo (section |5|). The only requirement on the hardware for our 
method to be efficient, is that the computer is quick in calculating elementary 
functions such as exp(a;), tan(x), and their inverse functions. Many present 
day computers which are equipped with co-processors for floating point calcu- 
lations are suitable for our methods. Our method which keeps uniformly high 
acceptance rate for all values of a, is particularly suitable to parallel or pipeline 
vector processors. 

From the general consideration given in section g, the use of the 'cosh' dis- 
tribution and the Cauchy distribution which we proposed as approximate dis- 
tributions, will be efficient for generating random variables taking values in a 
finite interval (e.g. [— tt, tt]) and whose distribution fa{d) is an even function 
and takes maximum a,t 6 — 0, and behaves like fa{0) ~ const. — a9^ near ~ , 



10 



with a parameter a > that controls the sharpness of the peak. This is a 
common feature of weight functions for statistical mechanical systems with one 
component spins and link variables. 

We would like to mention a couple of methods which appears in the litera- 
ture, both of which are based on rejection methods, but use different approxi- 
mate distributions / . The approximate distribution adopted in [^ is 

/>)=iV,exp(a(l--|0|)), 

TT 

where Na is a normalization constant. The acceptance rate R{a) for this choice 
is, from eq. (g|), 

1 aexp(— ca) 



R{a) 



Na exp(a) TT (1 — exp(— 2a)) ' 



2 2/2 

where Na is as in eq. (0) and c = — sin^^( — ) + Wl — ( — )^ — 1 « 0.2105 is a 

positive constant. In particular, at 'low temperature' a 3> 1 , 

2a 
R[a) ~ \ — exp(— c a) , 

which rapidly approaches as a —^ oo . Therefore this choice suffers from the 
same problem as with 'direct' method that as a ^ oo , the CPU-time to keep 
90% effective acceptance rate tends to oo . In fact, we have performed an effi- 
ciency test as explained in section ra for this choice of approximate distribution, 
and found that for a > 3 , the 'proposed cosh' method is faster. 

The approximate distribution of H is the Gaussian distribution which, in 

our notation, is 

/>)=7V,exp(-a(a)02)^ 



where 



1 /"^ 

-^ = exp(-a(a)6|2)d6l, 



2 1 

is a normalization constant, and a{a) — -^maxia,-). (To be precise, they 

TT"^ 4 

adopt 'direct' method for small a and Gaussian distribution for large a . We 
focus our attension on the large a case where the simple 'direct' method is not 
effective.) quotes Q for an algorithm of generating the Gaussian random 
variable. 

The original distribution fa{d) is now considered as a distributionon G R , 
where /a (6*) = if 16*1 > tt . The formula similar to eq. (g) (with — tt < 9 < tt 
replaced by € R) holds in the present case, and the acceptance rate R{a) for 
the Gaussian distribution is, 

1 
R{a) = 



Na exp(a) V TT^ 
11 



where Na is as in eq. (0). This coincides with the resuhs in H]. At a — > oo , 
R{a) — > — w 0.6366 . Since R{a) does not tend to zero, the effective CPU-time 
consumption is bounded for all range of a, which is of course a nice feature 
shared with our methods. The acceptance rate of 'proposed cosh' method is 
considerably higher than that of the Gaussian method for all values of a . As for 
a quantitative comparison of CPU-time consumption, H] compares the CPU- 
time between the 'direct' method and the Gaussian method, and finds that 
(with their machine) for a > 1.5 the Gaussian method becomes faster than the 
'direct' method, and slows down somewhat as a is increased further. (Roughly 
12% decrease in speed from a = 1.5 to a = cx) .) As can be seen from Figure 2, 
at a = 1.5 the 'proposed cosh' method is already faster (the ratio in speed is 
roughly 1.6) than the 'direct' method, and it practically does not slow down if 
a is increased. (A rough estimate shows that only 3% decrease in speed occurs 
from a = 1.5 to a = oo .) Therefore we may conclude that the 'proposed cosh' 
method is faster than the Gaussian method for a > 1.5 . 
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In this Appendix, we prove the statements in section 0. 

Let p be the density function for the distribution of the random variable (R): 

p{9) de = Prob[ h{Cji) (^[6,e + dO] ] , 

where 

Prob[P{i)]= lim —^{l<i<N\P{i)}. 

The statement p ~ fa is proved as follows. By eq. m) we have 

h-\0 + d0) - h-\0) + ^^^^2__d0 - h-\e) + md0, 

hence 

p{0) d0 = Prob[Cj, e [h-He),h-^0) + f{e) de] ] . (27) 
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Fix < a; < 1 . By definition we have < g{x) < 1 . Since oj'j , j = 1, 2, 3, • • • , 
is a uniform random variable, we see that the conditional acceptance rate for tui 
under the condition that LUi € [x,x + dx] , is g{x) . Therefore, 

ProbluJi e [x,x + dx] ] = Prob[LUi G [x,x + dx] ] x g{x) oc g{x) dx , 

where we also used the fact that w; , i = 1, 2, 3, • • • , is a uniform random variable. 
Insert this equation in eq. ( pTJ ) and use eq. (^ to obtain 

p{d) dd ex g{h-\e)) f{9) do ex fa{e) de . 

Since fa is correctly normalized, we have p{9) = fa{0)- This completes the 
proof. 

The proof that the acceptance rate is equal to R in cq. (pi) is as follows. 
For each x with < a^ < 1 , the conditional acceptance rate for tOi under the 
condition that u>i G [x,x + dx] , is g{x) . Therefore 



The acceptance rate ~ I g{x) dx 



Ri'm^dx 



/o fihix)) 
= R f k^f{e)d0 = R. 

J-TT f(0) 



This completes the proof. 



B 



In this Appendix, we give a proof that the 'optimal cosh' method in section 
correctly gives the random C/(l)-variables, and also we give the argument for 
the choice of the parameters. 

We consider the case a > a*. The proof for the case < a < a* is similar. 

By explicit calculation, one sees that ha^f^ of cq. ( p^ satisfies eq. (^. There- 
fore it suffices to show that eq. (jia) satisfies eq. (H) with eq. m). 

Define, 

0(9) = ^log<|4^^^ (28) 

^(cosh(a0) + /?)) . (29) 

(G depends on three free parameters a, a, and /?. We suppress the parameter 
dependences for the moment.) Then eq. (o) and eq. (m imply 

5(x) =exp{-a(G(/i(a;))- min G{e))} , (30) 

ee\-TT,Tr] 
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and 



R{a) = — ■^'■"^, , exp{a min G(e)} 



Comparing eq. pa) with eq. dlq), one sees that the results in section ra are 
correct if 

min G(0)=O. (31) 

ee[-7r,7r] 

Proposition. Let a > 0, a > 0, and /3 > —1. If G satisfies the three 
conditions 



then G satisfies 



G(a;)>0, -7r<a;<7r. (35) 



Assume for the moment that this Proposition is true. It is easy to see by 
explicit calculations that a = a(a) and /3 = /3(a) defined by eq. dl4) or eq. (la) 
satisfy the conditions (p2), (p3|), and (p4), and «(«) > and /3(a) > — 1 , for all 
a > a* . Since 

G(0) = 0, 

the Proposition implies that eq. ( pl| ) is satisfied for all a > a* . 
It remains to prove the Proposition. 

Proof of the Proposition. Since G is an even function, it is sufficient to 
prove G(x)> for < a; < tt . 

From (||) and (H) it follows that 0^(2 - /3) > /? + 1 . The equality holds 
if and only if a(/3 +1) ^ c? and o? = 3a — 1 , which, with ( |34| ) implies (3a — 
l)(exp(2a) — 1) > a(cosh(7rV3a — 1) — 1) . This is equivalent to a > a° (> 3/2) , 
where a° is defined by eq. (jl^). Therefore, if we define a set D by 

B = D1UD2, 

Di = {(a,a,/3) e (0,cx))^ X (-l,c»)| a^(2-/3) > /3 + 1 > aVa}, 

L>2 = {(a,a,/3) e (0,oo)2 X (-1,00)1 a^ =3a-l, /3+1 = aVa, a> 3/2}, 

it is sufficent to prove that for all (a, a, 0) d D and < x < tt , (B4) implies 
G(x) > . 
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Step 1. Fix (a, a, 0) £ D . Put g{x) = G'{x) . (This g has nothing to do with 



g in eq. (30).) Then we have 



fix) EE g{x)+g"(x) (36) 

a^sinh(ax) , / -2/ i / \ i\\ 

= —FT, TT — 'F^hia (cosh(Q:a;) - 1)) , 

a(/3 + cosh(aa;))3 ^ ^ ' " 

where 

h{y) = -y2 + (^ _ 2a-2(/3 + i))j^ + or\l5 + l)(a2(2 - /3) - (/3 + 1)) . 

Since (a, a, /3) G _D , we see that there exists one and only one positive root 
y = yo oi h{y) — and that 

h{y) > , if < y < yo , 
h{y) < , if y>yQ. 

Therefore if we let a; = xq to be the unique positive solution to the equation 
a~^(cosh(aa::) — 1) = yo , we have 

fix) > , if < a; < a;o , (37) 

fix) < , if x>xo. (38) 

Note that 5(0) = and ^'(0) > if (a, a,(3) eD. From eq. @) and eq. (|3) 
we therefore conclude 

5(2;) > , if < a; < xo and < a; < tt . (39) 



(The conclusion may be easily understood if one notes that eq. (36) is an equa- 
tion of motion of harmonic oscillation with external force / .) 
The equations (ll), (|3|), (|3|), imply 



G'{x) = g{x) > , if < a; < a;o and < a; < TT . (40) 

g{x) + g" (x) < , if a; > a;o . (41) 

Step 2. Fix a > and t = aa^^(/3+ 1) > 1 , and let a vary with the restriction 
(a, a, (3) e -D . The allowed region of a differs by the values of t and a: 

1. t > 1 , or f = 1 and a < ^7/2 . 

In this case, (a, a,l3) £ D is equivalent to a > (a^ + l)/3 . 

2. i = 1 and a > ^/7j2 . 

In this case, (a, a, /3) G -D is equivalent to a > {a^ + l)/3 . 
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Note that g{x) — ga{x) is continuous (uniformly continuous on compact sets 
in (0, tt] w.r.t. x) and increasing in a , and xo = xq^o) is continuous in a . Also, 

lini g{x) — sin a; 

a — *C30 

uniformly on compact sets in (0, tt] . 

We claim that for every a (such that (o, a, (3) ^ D), and for any xi and 3:3 
satisfying ga{xi) > 0, gaix^) > 0, and < xi < X3 < tt , we have ga{x) > 
0, X € [a;i,a;3] . 

Assume this is wrong: Assume that for a = uq and < xi < 0:2 < X3 < tt 
we have gao{xi) > M , gao{x2) < 0, and Qaoixa) > M , where M is a positive 
constant. Since gaix) is increasing in a , we have 

gaixi) > M , gaixa) > M , a > oq . 

Put 

q{a)= min ga{x) . 

Xi <X<X3 

Then q{a) is continuous in a and lim q{a) > . Therefore there exists ai > ao 

a — 'oo 

such that q{ai) = , which further implies that gai{x) > , xi < x < x^ ^ and 
that there exists X4 satisfying xi < 2:4 < x^ and gai{xi) = . In particular, 
gai{x4) — and g'^ (x^) > hold, which contradicts eq. ( [40|) and eq. ([4l|). 
Hence the claim is proved. 

Step 3. Fix (a, a, /?) G D . From the claim and eq. (Eol) we see that either 
g{x) — G'{x) > for < X < TT , or there exists x' such that < x' < tt , and 
g{x) > for < X < x' , and g{x) < for x' < x < tt . Hence G{x) is either 
increasing in < x < tt or has just one peak and no valley. Since G{0) = and 
G(7r) > , we have G{x) > , < x < tt . This completes the proof. 

We now turn to the argument for the choice of the parameters. We want to 
choose the parameters so that the acceptance rate is as large as possible. As 
stated at the end of section H, it is better to have as flat g{x) as possible, hence 
a flat G{x) . Thus we require G"(0) = 0. We impose the condition that the 
minimum of R{a) is achieved a.t 9 = . Note that this is equivalent to assuming 
eq. (|35|). As a necessary condition, we have G*^''-'(0) > and G(7r) > . (By the 
Proposition, we know that these are sufficient to ensure eq. (p3).) As we want 
to have flat G{x) , it should be best to have either G^^\0) = or G(7r) = . If 
one draws a graph of these three conditions, in (a, /3)-plane, one easily sees that 
the choice given in the section g] is the one that we are looking for. 

c 

In this Appendix, we give a sample FORTRAN program for generating random 
(7(1) variables, using 'proposed cosh' method. When the subroutine UlRND 
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is called with the parameter a as the first argument, it returns a random U{1) 
variable as the second argument. 

Note that this sample program is different from the program used to test 
the efficiency disscussed in section ra, where we modified the program in favor 
of a pipelined vector processor. 



SUBROUTINE U1RND(A,G) 
C VARIABLES FOR U(l) RANDOM NUMBERS 

REAL PI,A,AS,EPS,DEL,ALP,BET,DAP,DL1,DL2,BT1,H,G,GH,H1 

PARAMETER (PI=SNGL (3 . 14159265358979D0) ) 

PARAMETER (AS=SNGL (0 . 798953686083986D0) ) 

PARAMETER (EPS=0.001, DL1=0.35, DL2=1.03) 
C VARIABLES FOR UNIFORM RANDOM NUMBERS 

INTEGER C30 , CRND , IRND 

REAL RND 

REAL*8 C31 

PARAMETER (C30=2**30, C31=2D0**31, CRND=48828125) 

DATA IRND/1000001/ 
C 

DAP=MAX(0,A-AS) 

DEL=DLl*DAP+DL2*SqRT(DAP) 

ALP=MIN(SQRT(A*(2-EPS)), MAX(SqRT(EPS*A) ,DEL) ) 

BET=MAX(ALP*ALP/A, (C0SH(PI*ALP)-1)/(EXP(2*A)-1) )-l 

BT1=SQRT((1+BET)/(1-BET)) 
C 

1 CONTINUE 
C 

R=IRND/C31 

IRND=IRND*CRND 

IF (IRND .LT. 0) IRND=(IRND+C30)+C30 
C 

Hl=BTl*TAN((2*R-l)*ATAN(TANH(PI*ALP/2)/BTl)) 

H=ALOG ( ( 1+Hl) / (1-Hl) ) /ALP 

G=EXP (-A* ( 1-COS (H) ) ) * (COSH (ALP*H) +BET) / ( 1+BET) 
C 

R=IRND/C31 

IRND=IRND*CRND 

IF (IRND .LT. 0) IRND=(IRND+C30)+C30 
C 

IF (G .LT. R) GOTO 1 

RETURN 

END 
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Figure captions. 

Figure 1: Acceptance rate R{a); 'direct' method (lowest curve), 'optimized cosh' 
method (highest curve), and 'proposed cosh' mcthod(curve in between). 

Figure 2: Time consumption T to keep 90% acceptance rate R{a) > 0.9 . Iterated 
'direct' method (+) and 'proposed cosh' method (o). T is in aribtrary unit. 
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